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Adjusting the melting point of a model system via Gibbs-Duhem integration: 

application to a model of Aluminum 
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Model interaction potentials for real materials are generally optimized with respect to only those 
experimental properties that are easily evaluated as mechanical averages (e.g., elastic constants 
(at T=0 K), static lattice energies and liquid structure). For such potentials, agreement with 
experiment for the non-mechanical properties, such as the melting point, is not guaranteed and 
such values can deviate significantly from experiment. We present a method for re-parameterizing 
any model interaction potential of a real material to adjust its melting temperature to a value that 
is closer to its experimental melting temperature. This is done without significantly affecting the 
mechanical properties for which the potential was modeled. This method is an application of Gibbs- 
Duhem integration [D. Kofke, Mol. Phys.78, 1331 (1993)]. As a test we apply the method to an 
embedded atom model of aluminum [J. Mei and J.W. Davenport, Phys. Rev. B 46, 21 (1992)] 
for which the melting temperature for the thermodynamic limit is 826.4 ± l.SK - somewhat below 
the experimental value of 933K. After re-parameterization, the melting temperature of the modified 
potential is found to be 931. 5K ± 1.5K. 



I. INTRODUCTION 



The ability of a simulation to successfully predict the 
properties of real materials is primarily dependent upon 
the accuracy of the model interaction potential used. The 
construction of model interactions generally involves the 
optimization of the parameters of the potential with re- 
spect to mechanical properties of the material (crystal 
lattice constants, elastic constants, liquid density, etc.) 
as determined from experiment or ah initio calculations. 
Non-mechanical properties (i.e. those not obtainable as 
local averages over coordinates) such as phase transition 
temperatures are difficult to include in such optimization 
procedures and are generally calculated for the optimized 
model a posteriori, and the agreement of such quantities 
with experiment is not guaranteed. However, for some 
applications that deal directly with such properties, such 
as in the study of solid-liquid interfaces |2| in which the 
melting temperature plays an obviously important role, 
it is desirable to develop efficient procedures for including 
such non-mechanical properties in the optimization. In 
this work, we outline a general procedure for adjusting 
the potential parameters for a system designed to model 
a real system to improve the agreement of the melting 
point of that system with the experimental value. As an 
example we present an application to an existing embed- 
ded atom model of aluminum || . 

Our re-parameterization scheme is based on the power- 
ful Gibbs-Duhem integration method developed by Kofke 
QH. In this technique, the derivative along the coexis- 
tence curve of any coexistence condition (such as melting 



temperature or pressure) with respect to any parameter 
of the potential can be determined by an appropriate 
configurational average at a previously determined melt- 
ing point. The method is generated by the integration 
of a generalized Gibbs-Duhem equation and the steps 
are analogous to the derivation of the familiar Clapey- 
ron equation for the slope of the P-T coexistence curve. 
Gibbs-Duhem integration has been shown to be quite 
successful in efficiently determining the coexistence con- 
ditions for entire classes of potentials. For example, the 
phase diagram for the class of repulsive inverse power 
potentials, u{r) — e (— )™, was determined Q by starting 
with the known hard-sphere (n — oo) coexistence and in- 
tegrating the derivative of coexistence curve with respect 
to the parameter s = 1/n. The method has also proved 
useful in a variety of other applications 

In the current application, one begins with a model 
potential, parameterized for a real system in the usual 
way with respect to mechanical properties of the real 
system. The melting temperature (or pressure) for the 
model system is then calculated by thermodynamic inte- 
gration. Once this is done, the derivative of the melting 
temperature (or pressure) with respect to all parameters 
of the system can be determined via separate simula- 
tions on the coexisting fluid and solid using the Gibbs- 
Duhem procedure. The calculated derivatives allow an 
assessment of the effect of each individual parameter on 
the melting point. From this information an appropriate 
scheme to adjust the parameters to improve the melting 
point can be devised in such a way that the agreement 
with the other experimental properties is not unaccept- 
ably compromised. The Gibbs-Duhem integration and 
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our re-parameterization scheme is outlined in more de- 
tail in the next section. 

As a test application of this procedure, we examine an 
embedded atom model of aluminum developed by Mei 
and Davenport Q . This particular model was chosen for 
three reasons: First, the importance of aluminum as a 
material makes the development of an accurate model 
potential for simulation purposes desirable. Second, the 
large number of parameters and complicated nature of 
the embedded atom potential increases the need for a 
systematic, as opposed to ad hoc procedure for adjusting 
the melting point. In addition, the zero-pressure melting 
point for the Mei-Davenport potential has been previ- 
ously determined 0] to be 800±9K somewhat lower than 
the experimental melting point of aluminum at 933K. 
(Note that, the melting point determined by Mei and 
Davenport was calculated for a 256 particle system - the 
actual value for this potential in the thermodynamic limit 
is slightly higher at 826.4±1.3K.) We find that, in this 
specific case, only one of the parameters of the poten- 
tial has any significant effect on the potential and that 
changing this parameter according to the Gibbs-Duhcm 
procedure yields a new model with the correct experi- 
mental melting point with no significant change in the 
quality of the agreement of the quantities with respect 
to which the model was originally optimized. Details of 
this calculation as well as a description of the model can 
be found in Section || below. 



II. GIBBS-DUHEM INTEGRATION AND 
MELTING POINT OPTIMIZATION 

The technique of Gibbs-Duhem integration has been 
well described previously by Kofke ||| , but in the interests 
of completeness and clarity we repeat the basic derivation 
here. Consider a single-component system with an arbi- 
trary interaction potential, U({Ri}, {Xi}), where the Ri 
are the atomic coordinates and the Xi are the parameters 
that define the potential - no restriction to pairwise ad- 
ditivity need be assumed. Assume there are two phases 
a and (3 in coexistence at a temperature T and pressure 
P. On the surface of coexistence, the chemical poten- 
tials (molar Gibbs free energies) of the two phases must 
be equal. To quantify how changes in P,T, and X^ will 
affect the chemical potential one can define a generalized 
Gibbs-Duhcm equation, 



-sdT + vdP 



(1) 



where \i is the chemical potential, s and v are the entropy 
and volume per particle, respectively, and the Xi are gen- 
eralized thermodynamic variables conjugate to the poten- 
tial parameters, Xi, defined as 



NX % 



f 8G 

\dXi j T,p,x jtj ^i 



(2) 



Now as one moves infinitesimally away from the origi- 
nal coexistence point (P,T, {Xi}) to another point (P + 
dP, T + dT, {Xi + dXi}) on the surface of coexistence, the 
change in fi must be identical in both phases. This con- 
dition together with Eq. [j] gives 

- fi/3 = "(Sq - Sp)dT + (v a - Vp)dP 

+ ^( A «.i-A/3,i)c^Q = 0, (3) 

i 

where fi a and fj,p are the chemical potentials for each 
of the respective phases. Assuming constant pressure 
(dP = 0), since we are interested here in changes in the 
transition temperature, the previous equation can be re- 
arranged to give 



dT 
dX~ 



T(X a i — Xf3 z ) _ TAXi 



hp 



Ah 



(4) 



where we have also assumed that at coexistence, As = 
Ah j T, where Ah is the latent heat per particle for the 
phase transition. (Note that, the corresponding equa- 
tions for the ( j^- J can be easily obtained by 

replacing Ah/T in Eq. Ejwith Aw.) 

The Xi can be related to mechanical averages that can 
be easily calculated in a molecular dynamics or Monte 
Carlo simulation. First, the Gibbs free energy is related 
to the isothermal-isobaric distribution, A(N,P,T) , as 
follows. 



G = -k b T\nA(N, P, T) 



(5) 



which for a classical system with interaction potential 
U{{Yi}) is given by 

A{N ' P ' T) = A^AM dV J CXP (_/? U ~ PPV) ' 

(6) 

where V is the volume. Taking the derivative of Eq. [| 
with respect to the parameter Xi gives 



dG 
dX~ 



= -kT 



<91nA 



dX t 



* / T,P,X, 



1 



- / dV I d N r 



8U 

io J \')X./ x :/ .., 

x exp (-pU - (3PV) . (7) 



Using Eq. we have 
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' dG 

dXr 



T,P.X, 




(8) 



N,P.T 



Using Eq. |5] and Eq. || we can calculate the change 
necessary in each of the parameters {A^}, to alter the 
melting temperature of our interaction potential by some 
arbitrary amount. As long as the changes are not too 
large, the calculation can be greatly simplified by the 
assumption that the deviation remains linear. 



dT 
dA~ 



(Xi — Xifi) 



Xa , jj,coex 



(9) 



where T m ^ is the melting point of the original model with 
parameters {A^o}- (If the linear approximation does not 
hold - and it should always be checked - one could inte- 
grate the differential equation using an appropriate nu- 
merical technique.) From a single simulation, one could 
calculate the necessary change in all the different parame- 
ters of the potential corresponding to a particular change 
in the melting temperature. One could then, use this in- 
formation to construct a cost function to re-optimize the 
potential that includes the previously used set {A^ x ' 3 '} 
of experimentally determined (mechanical) properties as 
well as the experimental melting temperature T^^' as 
optimization targets. For example, a linear least squares 
procedure could be utilized using appropriately chosen 
set of weight functions Wf. 



,/({X 4 }) = u; T [^ x P--T m ({Xj)] 2 

+ E w 'H exp ' " Mix,})? 



(10) 



(In the example application discussed in the next section, 
we found that only one of the several parameters of the 
potential had any significant effect on the melting point, 
which greatly simplified the re-parameterization by mak- 
ing it possible to modify the potential without having 
to directly use such a cost function. It is not expected, 
however, that this will be true, in general.) 



III. RESULTS FOR A MODEL OF ALUMINUM: 
DETERMINING THE ORIGINAL MELTING 
POINT 

As an application of the method we examine an em- 
bedded atom model |l^] for aluminum developed by Mei 
and Davenport ||. The parameters in their potential 
were fit in order to optimize the potential with respect 
to a variety of experimental properties such as the co- 
hesive energy (E c ), lattice constant (v^ro), un-relaxed 



vacancy-formation energy, and elastic constants of the 
static fee crystal at zero temperature. The total poten- 
tial energy has the form 



(11) 



where 

F(p) = -E c 



*/0 
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-^00 E s mexp[-(Vm- 1)7] 
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Pr = E /0 
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(r) - -<h 
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exp 
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(12) 

(13) 

(14) 
(15) 



and is used in conjunction with the following cutoff func- 
tion: 

1 r < r n , 

q(r) = { (1 - x) 3 (l + 3x + 6x 2 ) r n < r < r c , (16) 

r > r c 



x = (r - r n )/(r c - r n ) 



(17) 



with a cutoff distance (r c ) between the third and fourth 
neighbor shells of a static fee crystal. Using values of 
r„ = 1.75ro and r c = 1.95ro , the functions q(r)f(r) and 
q(r)4>(r) go smoothly to zero at r = r c . The parameters 
of the potential fp|]Tl]] are given in Table |. For this em- 
bedded atom model, mass is measured in amu, distance 
in A, and energy in eV. The natural simulation time unit 
is calculated to be 10.181fs. 

To begin, the melting temperature of the original 
model (at P = 0) must be determined. Mei and Dav- 
enport perform a calculation of the free energies of the 
aluminum melt and fee crystal using thermodynamic in- 
tegration of the Gibbs free energy. For a 256 particle 
system, they obtain a melting temperature of 800±9K. 
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Table |: Parameters for the Aluminum embedded atom potential developed by Mei and Davenport. 





E c 




ro 


a 


P 


7 


<5 


3.39 


0.1318 


2.8638 


4.60 


7.10 


7.34759 


7.35 


CO 


Cl 


C2 




C4 


C5 




0.64085 


-6.83764 


26.75616 


-47.16495 


36.18925 


-8.60834 







TABLE I. E c and <^> are in units of eV. r is in units of A, and the other parameters are dimensionless. 



In repeating their melting point determination, we 
found that the value they obtained is not quite correct, 
due primarily to the small system size studied and a prob- 
lem in the choice of thermodynamic integration path for 
the liquid phase. Our melting point determination was 
performed using the same basic methodology as Mei and 
Davenport, described below. 

The Gibbs free energies per particle of the liquid and 
crystal as a function of temperature at constant pressure 
can be obtained by thermodynamic integration using the 
integral form of the Gibbs-Hclmholtz equation 



g(T,P) g(T ,P) 



T 



To 



e(r,P) 



dr , 



(18) 



where T is a predetermined reference temperature and 
e(T, P) is the average total energy per particle. The 
Gibbs free energy at the reference temperature must be 
obtained separately by thermodynamic integration from 
a suitable ideal reference state. To do this, the interaction 
potential is parameterized along a linear path between 
that of the reference potential, Uq and the full potential 
U 



W(0=#/+(1-O«b 



(19) 



The Gibbs free energy per particle relative to that of the 
reference system can then be obtained by thermodynamic 
integration along the path: 

g(T,P) = g(T,P;Z=l) 

= g(T ,P^ = 0) + J\^^-^. (20) 

For the crystal, the reference system chosen is that of 
an Einstein crystal |L2[ 



U E {{n}) = \mu? D 5^(r< - r„) 2 , 



(21) 



where the {r^o} are the ideal crystal lattice positions, m 
is the mass and u>n is the Debye frequency, which to mini- 
mize the difference between the reference and full system, 



should be chosen to give a similar mean squared displace- 
ment for the atoms at the temperature of interest. For 
this system at 296K, the optimum Debye frequency corre- 
sponds to a Debye temperature (To = huo/k) of 207K. 
Note that this frequency is different than the one used 
by Mei and Davenport. 

The reference system for the liquid phase is an ideal 
gas, but the transformation must performed as a two 
step process in order to avoid the liquid/gas phase transi- 
tion. Mei and Davenport follow the reversible expansion 
method used by Broughton and Li [jl3|. The two step 
process is accomplished by first turning off the attrac- 
tive part of the potential followed a volume expansion 
to reach the ideal gas limit. In step one, it is extremely 
important to turn off the attractive part of the poten- 
tial in a way that will not drastically alter the effective 
radius of the potential. If care is not taken, the system 
will freeze during the integration of this step. Here, we 
write the interatomic potential as a linear interpolation 
between the actual potential and the reference system 



<f>(r; C) = <j)re P {r) + ( <j) a tt(r) 



(22) 



As the value of £ is varied from 1 to 0, the system is 
transformed from the original Mei and Davenport poten- 
tial to that of a purely repulsive potential. Step two of 
the integration is a volume expansion. The free energy 
change in this step is given by 



A t T f d P' 

Jo P 



(3P 
~7 



- 1 



(23) 



The potential splitting that Mei and Davenport use has 
the problem in that their repulsive potential 4> rep = 



o<5 exp -7 (j^ - l) 



has an effective radius that is 

much larger than the actual potential (See Fig. 0a) so 
that as the rest of the potential is turned off the sys- 
tem freezes as the effective packing fraction increases. To 
remedy this we use a WCA |lj] splitting where the re- 
pulsive part of the potential is equal to the the potential 
energy for radii less than the radius at the minimum of 
the potential and zero for larger radii. This prescription 
(illustrated in Fig. |l|) gives an effective radius more sim- 
ilar to the full system and avoids the freezing transition. 
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FIG. 1. The different splitting methods for the pair part 
of the aluminum embedded atom potential are graphed sep- 
arately as a function of r. (a) Mei and Davenport splitting 
(b) WCA splitting method. Note the difference in the effec- 
tive radii between the combined and repulsive parts of the 
potential. In the Mei and Davenport splitting this difference 
is large and can lead to freezing during the course of integra- 
tion, whereas the WCA method has a much smaller difference 
in effective radii. 

In order to obtain energy curves needed in Eq. [l8] 
as a function of temperature, molecular dynamics (MD) 
simulations were conducted at several different temper- 
atures using the Nose-Poincare- Anderson algorithm |l5| 
for isothermal-isobaric molecular dynamics. The fee crys- 



tal was simulated at 50K intervals from 296K to 946K, 
while the liquid was studied over a smaller temperature 
range from 762K to 1152K using 30K intervals. All of 
the isothermal-isobaric MD simulation runs were equili- 
brated for 100,000 steps and sampled for 300,000 steps. 
From these simulations we obtained the average energy 
and density for both the crystal and liquid as a func- 
tion of temperature and system size. Both the energies 
and densities were fit to second order polynomials. The 
coefficients for these polynomials are shown in Table EL 
These polynomials were used in the construction of the 
free energy curves as described in Eq. [j^. 
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FIG. 2. Integrand for the numerical integration for the 
solid free energy at 296K as a function of £. Values of £ were 
chose based on the ten-point Gauss-Legendre quadrature. As 
the value of £ changes from 1 to 0, the system is transformed 
from an embedded atom solid to that of an Einstein crystal. 

Using the density obtained from the constant NPT 
simulations, the Gibbs free energy of the fee crystal at 
To = 296AT was calculated by running several simula- 
tions at constant N,V,T using the Einstein crystal refer- 
ence state to evaluate the integrand of Eq. ^o] at values of 
£ corresponding to a ten-point Gauss-Legendre quadra- 
ture. The simulations were performed using the Nose- 
Poincare algorithm Figure || shows a plot of the 
integrand vs. £ for all of the system sizes studied. At a 
value of £ — 1 the system is governed solely by the 
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Table [il]: Coefficients for the Energy and Density Polynomials 



N 


256 




500 


864 




2048 


4000 




Liquid Energy ao 


-3.95773 x 10 


-8 


-2.62765 x 10~ s 


-3.47219 x 10 


-8 


-2.82724 x 10 -i< 


-3.14663 x 10 


-8 


«i 


4.24388 x 10" 


4 


3.98266 x 10- 4 


4.15658 x 10- 


4 


4.02841 x 10 -4 


4.08950 x 10- 


4 


0.2 


-3.39543 




-3.38194 


-3.39085 




-3.38448 


-3.38731 




Solid Energy ao 


5.75263 x 10" 


8 


5.85627 x 10" 8 


5.88951 x 10- 


8 


5.88624 x 10~ 8 


5.88849 x 10" 


8 


ai 


2.30701 x 10" 


4 


2.30331 x 10 -4 


2.30283 x 10- 


4 


2.30528 x 10" 4 


2.30565 x 10" 


4 


a 2 


-3.38550 




-3.38533 


-3.38529 




-3.38532 


-3.38531 




Liquid Density ao 


7.59175 x 10 - 


10 


1.61803 x 10- iu 


5.35909 x 10" 


10 


3.19601 x 10- 1U 


3.91853 x 10" 


1U 


Oi 


-1.03177 x 10 


-5 


-9.15213 x 10- 6 


-9.91989 x 10 


-6 


-9.47974 x 10~ 6 


-9.62320 x 10 


-6 


a 2 


6.09309 x 10 _ 


2 


6.03495 x 10" 2 


6.07444 x 10" 


2 


6.05232 x 10 -2 


6.05920 x 10" 


2 


Solid Density ao 


-2.59267 x 10 


-9 


-2.62088 x 10- 9 


-2.63749 x 10 


-9 


-2.64149 x 10~ 9 


-2.64086 x 10 


-9 


a x 


-2.79768 x 10 


-6 


-2.78739 x 10 -6 


-2.77752 x 10 


-6 


-2.77869 x 10" 6 


-2.78125 x 10 


-6 


a 2 


6.00427 x 10- 


2 


6.00404 x 1Q- 2 


6.00379 x 10- 


2 


6.00380 x 1Q- 2 


6.00387 x 10- 


12 



TABLE II. Coefficients for the solid and liquid energy and density curves. The curves are polynomials of the form: 
y — a^x 2 ± a\x ± a 2 . 



embedded atom potential. As the value of £ goes to zero, 
the potential is gradually changed to that of an Einstein 
crystal. For each value of £, the system was equilibrated 
for 100,000 steps and sampled for 100,000 steps. 

Simulations for step one of the liquid free-energy at a 
temperature of 1092.fr were performed in the same man- 
ner as for the crystal. Here the attractive part of the 
potential (Eq. 22) was slowly turned off as the value £ 
was changed from one to zero. This is shown graphi- 
cally in Figure ||(a). As with the solid, ten-point Gauss- 
Legendre integration was used to numerically compute 
this integral. The second step of the liquid free-energy 
calculation included a series of constant NVT simula- 
tions at decreasing densities starting with the repulsive 



potential system from the conclusion of step one. Aver- 
age values for the step two integrand (Eq. |2^) are shown 
graphically in Figure ||(b) . The simulations for both steps 
of the liquid free energy integrations were equilibrated 
for 100,000 steps and sampled for 100,000 steps at each 
value of the integrand. The numerical integration for step 
two was performed using the ten-point Simpsons quadra- 
ture. This method was chosen over the Gauss-Legendre 
quadrature due to the inaccuracy of the sampling at low 
density. The value of the integrand at zero density was 
obtained by an analytical calculation of the second virial 
coefficient (i?2)- 

Free energies and the melting point are shown in Ta- 
ble [II as a function of particle number. A graph of the 



Table [II: Free Energies and Melting Points for 8 = 7.35 



N 


256 


500 


864 


2048 


4000 


Original 5 


7.35 


7.35 


7.35 


7.35 


7.35 


9s{T 


= 296K) 


(eV/atom) 


-3.40200(3) 


-3.40219(2) 


-3.40236(2) 


-3.40246(1) 


-3.40254(1) 


9l(T = 


1092A") 


(eV/atom) 


-3.8561(6) 


-3.8561(4) 


-3.8563(3) 


-3.8565(2) 


-3.8565(2) 






T m (K) 


802.8±5.6 


814.3± 3.7 


819.7±2.9 


822.5±1.9 


825.2±1.3 


New 8 


8.70 


8.57 


8.50 


8.47 


8.45 


9s{T 


= 296if) 


(eV/atom) 


-3.39538(3) 


-3.39615(2) 


-3.39664(1) 


-3.39687(1) 


-3.39704(1) 


9l(T = 


1092iO 


(eV/atom) 


-3.8135(6) 


-3.8175(5) 


-3.8200(3) 


-3.8209(2) 


-3.8217(2) 






T m (K) 


934.9 ± 5.9 


933.1 ±4.4 


930.9 ± 3.0 


931.9 ±2.0 


931.5 ±1.5 



TABLE III. Simulation averages for the solid and liquid Gibbs free energies along with the calculated melting temperature 
for several system sizes. Free energies are in units of eV/atom and the melting temperature is in K. 



G 



melting temperature vs. 1/N (Figure [!]) shows that 
at infinite particle number the melting point for the em- 
bedded atom potential proposed by Mei and Davenport 
approaches 826. AK ± 1.3K. It should be noted that the 
major error in the melting point calculated by Mei and 
Davenport for their potential 800±9if is primarily due 
to the small size of their system. The problems with 
their potential splitting appears to have had little effect, 
probably due to cancellation of errors, as the melting 
point that we determine here for the 256 particle system 
agrees with theirs within the simulation error. Recently, 
Morris, Wang, Ho and Chan argued on the basis of 
the stability of crystal-liquid interfaces that the melting 
point of the Mei-Davenport aluminum potential is actu- 
ally significantly lower (around 725K). Our results do not 
support this conclusion and as a check we have carefully 
set up stable stress-free interfaces at our calculated melt- 
ing point. The systems exhibit melting (freezing) as the 
temperature is raised (lowered) away from our calculated 
melting point. The lower temperature transitions of Mor- 
ris, et al. were most likely due to significant un-relaxed 
stress in the crystal. 



IV. RESULTS FOR A MODEL OF ALUMINUM: 
RE-PARAMETERIZING THE POTENTIAL 

To adjust the melting temperature of the embedded 



to a melting temperature of T m — 933i-T. These values 



are listed in Table [II. (Note that, although we did an 
accurate numerical integration along the coexistence 



for 



atom potential for aluminum, we first calculate -^-r 
each of the parameters in the potential. (Note, for sim- 
plicity and consistency the expansion coefficients, q, in 
Eq. |lj were kept constant). The derivatives were cal- 
culated in a single (constant NPT) simulation at the 
melting temperature calculated in the previous section 
with P = 0, for each of several system sizes (N =256, 
500, 864, 2048 and 4000). The system was equilibrated 
for 100,000 steps followed by 100,000 steps for average- 
ing. From this simulation it was determined that 
was significant only for the parameter 8 - the other pa- 
rameters of the potential have little effect on the melting 
point. The complicated nature of the potential makes 
it difficult to assign any physical explanation to the sen- 
sitivity of the melting point to 6 relative to the other 
parameters. 

Next, a series of simulations were performed for each 
system size to integrate along the coexistence curve from 
the initial calculated melting point of the potential to the 
true experimental melting temperature. At each temper- 
ature along the coexistence curve, the system was equili- 
brated for 100,000 steps and sampled over 300,000 steps. 
In our experiments we use a fourth-order predictor cor- 
rector integrator to carry out the integration along the 
coexistence curve as a function of S. Fig shows this in- 
tegration graphically with final values of d corresponding 




-3.200 



FIG. 3. (a) Simulation results for the integrand of step 
one of the liquid free energy calculation at T = 1092_K" as 
a function of £. This integration slowly turns off the at- 
tractive part of the potential as the value of £ changes from 
1 to 0. Values for £ were again based upon the ten-point 
Gauss- Legendre quadrature, (b) Integrand for the volume 
expansion integration (step 2) in the liquid free energy calcu- 
lation at T = 10927f as a function of p* — — . The numerical 

r PO 

integration for this step was performed using the ten-point 
Simpson quadrature. The value at p* — was obtained by an 
analytic calculation of the second virial coefficient (i?2(T)). 
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curve, the results indicate that using the approximation 
that the deriviative is a constant in the region of interest 
would have been correct to within the simulation error.) 



Table [II. The experimental value for the melting tem- 



830 




0.001 



0.002 



1 / N (atoms 



0.003 



0.004 



FIG. 4. The melting temperature of the embedded atom 
potential is plotted as a function of inverse particle size. As 
N goes to infinity (i —> 0), the melting temperature ap- 
proaches 826. 4K. 

In order to confirm that the melting point did indeed 
change as expected, the melting point calculation was re- 
peated using the newly calculated value of S. The new 
melting temperature for the embedded atom potential 
with a value of <5 = 8.45 corresponding to a system size 
of N=4000, was calculated to be 931.5^ ± 1.5K. Melt- 
ing temperatures for the other system sizes are listed in 



perature of aluminum is 933. 47K. 

Mei and Davenport initially determined 6 (and the 
other parameters in the potential) by fitting the poten- 
tial to certain physical constants. For the new value of S 
(8.45), we have re-calculated a variety of physical prop- 
erties of aluminum. These quantities, for the original 
potential, the reparameterized potential and the experi- 
mental values are collected in Table [V. From this data, 




«J 8.00 - 



7.75 - 



7.50 



7.25 



800 



820 



S40 



860 



880 



900 



920 



940 



Melting Temperature (K) 



FIG. 5. Values of 5 along the Gibbs-Duhem integration of 
the crystal-melt coexistence curve. These curves were gen- 
erated using a fourth-order predictor-corrector and show the 
melting temperature as a function of 8 for each of the various 
system sizes. 



Table IV: Physical constants for N=4000 



5 = 7.35 8 = 8.45 Experimental 



Un-relaxed Vacancy Formation (eV/atom) 
Latent Heat (eV/atom) 
Cii (dynes/cm 2 ) 
C12 (dynes/cm 2 ) 
C44 (dynes/cm 2 ) 



4.07211 


4.18679 


N/A 


0.0830 


0.0973 


0.111 


0.93 


0.96 


1.07 


0.69 


0.68 


0.61 


0.33 


0.36 


0.28 



TABLE IV. Here the vacancy formation energy, latent heat values and elastic constants are presented for the original version 
of the EAM Al potential and the modified version. C\i and C12 for the modified version (with the experimental melting 
temperature) are closer to their experimental values. C44 seems to get slightly worse. 



8 



we see that, in comparison to the original potential with 
S = 7.35, the new potential more closely models the ex- 
perimental values of the (T = OK) elastic constants 
Cii and C22 while C44 becomes slightly worse in compar- 
ison to its experimentally determined value. In addition, 
mostly as a consequence of the improved melting point, 
the latent heat is considerably improved. 



V. SUMMARY 

We have outlined an application of the Gibbs-Duhem 
integration method of Kofke Q, with which a model 
interaction potential can be re-parameterized, including 
the experimental melting point in the optimization proto- 
col. The melting temperature of a potential can then be 
adjusted similar to the tuning of other parameters in the 
potential. Since non- mechanical properties, such as the 
melting point, are not generally included in potential op- 
timization and the agreement of such quantities with ex- 
periment is not guaranteed, such a procedure will be use- 
ful in situations, such as in the simulation of crystal-melt 
interfaces, where the having the correct melting point is 
highly desirable. The method is general and can easily be 
extended to a variety of systems. As an example of the 
utility of the method, we apply the procedure to repa- 
rameterizing a popular model of aluminum Q for which 
the melting point has been calculated to be over 100K 
below the accepted experimental value. We demonstrate 
that the reparamterized potential has a melting point 
that agrees within the statistical error with the experi- 
mental value of 933K and that reparameterization does 
not degrade the quality of the potential with respect a 
variety of properties, and that, in fact, quantities such 
as the elastic constants Cn and C12 and the latent heat, 
agreement is improved. 
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